Figure 2: Maps of Species CHI

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, fig.height = 4, fig.width = 7)

library(terra)
library(sf)
library(oharac)
library(data.table)
library(tidyverse)
library(biscale)  ### for bivariate legend
library(cowplot)
library(here)
source(here('common_fxns.R'))

1 Summary

Create Fig. 2 for manuscript: A) a panel of cumulative impacts, B) & C) separated into climate and non-climate, D) quartile overlaps

2 Data

Results from scripts in folder 2_process_spp_vuln_and_impacts

3 Methods

Read in the output rasters for species-level impacts (unweighted) - total, climate, and non-climate; reclassify based on quartile/quintile for biplot.

ocean_map <- rast(here('_spatial/ocean_area_mol.tif'))
ocean_df <- as.data.frame(ocean_map, xy = TRUE)

land_sf <- read_sf(here('_spatial/ne_10m_land/ne_10m_land_no_casp.shp')) %>%
  st_transform(st_crs(ocean_map))

nspp_mask <- rast(here('_output/nspp_maps/species_richness.tif'))
nspp_vals <- values(nspp_mask)
values(nspp_mask)[nspp_vals < 20] <- NA

chi_gps_spp <- list.files(here('_output/cumulative_impact_maps'),
                         pattern = 'species_group_chi',
                         full.names = TRUE)

climate_map_spp     <- rast(chi_gps_spp[1])

nonclimate_map_spp  <- rast(chi_gps_spp[2:4]) %>%
  sum(na.rm = TRUE)

chi_total_map_spp <- rast(here('_output/cumulative_impact_maps/chi_species.tif'))

3.1 Create panels

3.1.1 Panel A: total CHI

Squashing extreme outliers!

squash_outliers <- quantile(values(chi_total_map_spp), .9999, na.rm = TRUE)

cc_spp_df <- as.data.frame(climate_map_spp, xy = TRUE) %>%
  setNames(c('x', 'y', 'cc_spp')) %>%
  mutate(cc_spp = ifelse(cc_spp < squash_outliers, cc_spp, squash_outliers))
noncc_spp_df <- as.data.frame(nonclimate_map_spp, xy = TRUE) %>%
  setNames(c('x', 'y', 'noncc_spp')) %>%
  mutate(cc_spp = ifelse(noncc_spp < squash_outliers, noncc_spp, squash_outliers))
chi_spp_df <- as.data.frame(chi_total_map_spp, xy = TRUE) %>%
  setNames(c('x', 'y', 'chi')) %>%
  mutate(chi = ifelse(chi < squash_outliers, chi, squash_outliers))

chi_lims <- c(0, round(max(chi_spp_df$chi), 2))
chi_brks <- c(0, .25, .5, 0.75, round(max(chi_spp_df$chi), 2))
chi_lbls <- sprintf('%.2f', chi_brks)
chi_lbls[length(chi_lbls)] <- paste0('≥ ', chi_lbls[length(chi_lbls)])

xfm <- function(x) {x^1} ### function to warp CHI scale for visual impact

globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                      c(180, 90), c(180, -90), c(-180, -90)) 
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = crs(ocean_map))

inset_bbox_df <- tribble(
  ~xmin,   ~xmax,  ~ymin,  ~ymax, ~loc,
 -1.0e6,  +3.5e6, +3.5e6, +8.0e6, 'europe',
 +8.5e6, +13.5e6, -0.0e6, +5.0e6, 'se_asia',
 -9.5e6, -5.5e6, +0.5e6, +4.5e6, 'caribbean'
 )

panel_a_chi <- ggplot() +
  ### plot data:
  geom_raster(data = chi_spp_df, aes(x, y, fill = xfm(chi))) +
  ### continents:
  geom_sf(data = land_sf,
          fill = 'grey96', color = NA, 
          size = .1) +
  geom_sf(data = globe_border,
          fill = NA, color = 'grey70', 
          size = .1) +
  geom_rect(data = inset_bbox_df, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            color = 'red', fill = NA, linewidth = 0.2) +
  scale_fill_viridis_c(limits = xfm(chi_lims), breaks = xfm(chi_brks), labels = chi_lbls) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_void() +
  theme(axis.text = element_blank(), 
        axis.title = element_blank(),
        legend.text = element_text(hjust = 0.5, vjust = 1, size = 7, angle = 0),
        legend.title = element_blank(),
        legend.key.height = unit(.16, 'in'),
        legend.key.width = unit(.32, 'in'),
        legend.direction = 'horizontal')

# panel_a_chi
# asdf <- get_legend(panel_a_chi)

plot_inset <- function(chi_df, bbox_df) {
  ### for chi_df, rename the chi column to just 'chi'
  ### for bbox_df, filter to the appropriate row
  ggplot() +
    ### plot data:
    geom_raster(data = chi_df, aes(x, y, fill = xfm(chi)), show.legend = FALSE) +
    scale_fill_viridis_c(limits = xfm(chi_lims), breaks = xfm(chi_brks), labels = chi_lbls) +
    scale_x_continuous(expand = c(0, 0), limits = c(bbox_df$xmin[1], bbox_df$xmax[1])) +
    scale_y_continuous(expand = c(0, 0), limits = c(bbox_df$ymin[1], bbox_df$ymax[1])) +
    theme_void() +
    theme(panel.border = element_rect(color = 'black', fill = NA, linewidth = .5),
          panel.background = element_rect(color = NA, fill = 'grey96'))
}

inset_a_i <- plot_inset(chi_df = chi_spp_df, bbox_df = inset_bbox_df %>% filter(loc == 'europe'))
inset_a_ii <- plot_inset(chi_df = chi_spp_df, bbox_df = inset_bbox_df %>% filter(loc == 'caribbean'))
inset_a_iii <- plot_inset(chi_df = chi_spp_df, bbox_df = inset_bbox_df %>% filter(loc == 'se_asia'))

panel_b_cc <- ggplot() +
  ### plot data:
  geom_raster(data = cc_spp_df, aes(x, y, fill = xfm(cc_spp))) +
  ### continents:
  geom_sf(data = land_sf,
          fill = 'grey96', color = NA, 
          size = .1) +
  geom_sf(data = globe_border,
          fill = NA, color = 'grey70', 
          size = .1) +
  # geom_rect(data = inset_bbox_df, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
  #           color = 'red', fill = NA, linewidth = 0.2) +
  scale_fill_viridis_c(limits = xfm(chi_lims), breaks = xfm(chi_brks), labels = chi_brks) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_void() +
  theme(axis.text = element_blank(), 
        axis.title = element_blank())

# panel_b_cc

panel_c_noncc <- ggplot() +
  ### plot data:
  geom_raster(data = noncc_spp_df, aes(x, y, fill = xfm(noncc_spp))) +
  ### continents:
  geom_sf(data = land_sf,
          fill = 'grey96', color = NA, 
          size = .1) +
  geom_sf(data = globe_border,
          fill = NA, color = 'grey70', 
          size = .1) +
  # geom_rect(data = inset_bbox_df, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
  #           color = 'red', fill = NA, linewidth = 0.2) +
  scale_fill_viridis_c(limits = xfm(chi_lims), breaks = xfm(chi_brks), labels = chi_brks) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_void() +
  theme(axis.text = element_blank(), 
        axis.title = element_blank())

# panel_c_noncc

3.1.2 Panel D: spp climate vs nonclimate bivariate plot

Use a four-level bivariate color scale to capture spatial overlap of quartiles for climate and non-climate stressors.

Reclassify values using ntile to get quartile values:

cc_qtile_spp <- noncc_qtile_spp <- climate_map_spp ### initialize two new maps
values(cc_qtile_spp) <- ntile(values(climate_map_spp), 4)
values(noncc_qtile_spp) <- ntile(values(nonclimate_map_spp), 4)

cc_spp_qtile_df <- as.data.frame(cc_qtile_spp, xy = TRUE) %>%
  setNames(c('x', 'y', 'cc_spp'))
noncc_spp_qtile_df <- as.data.frame(noncc_qtile_spp, xy = TRUE) %>%
  setNames(c('x', 'y', 'noncc_spp'))

cc_noncc_spp_qtile_df <- full_join(cc_spp_qtile_df, noncc_spp_qtile_df, by = c('x', 'y')) %>%
  mutate(bi_class = paste(cc_spp, noncc_spp, sep = '-'))

cc_noncc_summary <- cc_noncc_spp_qtile_df %>%
  group_by(cc_spp, noncc_spp) %>%
  summarize(n = n(),
            pct = n() / nrow(cc_noncc_spp_qtile_df), 
            .groups = 'drop') %>%
  mutate(pct = round(pct * 100, 1),
         lbl = paste0(pct, '%'))

knitr::kable(cc_noncc_summary, digits = 3)

cc_spp	noncc_spp	n	pct	lbl
1	1	324581	8.8	8.8%
1	2	300192	8.2	8.2%
1	3	179527	4.9	4.9%
1	4	115698	3.1	3.1%
2	1	279191	7.6	7.6%
2	2	250271	6.8	6.8%
2	3	231144	6.3	6.3%
2	4	159392	4.3	4.3%
3	1	211177	5.7	5.7%
3	2	178541	4.9	4.9%
3	3	264828	7.2	7.2%
3	4	265452	7.2	7.2%
4	1	105049	2.9	2.9%
4	2	190994	5.2	5.2%
4	3	244499	6.6	6.6%
4	4	379456	10.3	10.3%

For the bivariate color plot, one axis of colors represents climate, the other represents non-climate; for corner values (max agreement and max disagreement), use diverging color palettes; for intermediate values use a light shading to indicate variation more subtly.

map_pal_raw <- bi_pal(pal = 'PurpleOr', dim = 4, preview = FALSE)
  
map_pal_mtx      <- matrix(map_pal_raw, nrow = 4, ncol = 4)
map_pal_mtx[3, ] <- colorspace::lighten(map_pal_mtx[3, ], .1)
map_pal_mtx[2, ] <- colorspace::lighten(map_pal_mtx[2, ], .2)
map_pal_mtx[1, ] <- colorspace::lighten(map_pal_mtx[1, ], .3)
map_pal_mtx[ , 3] <- colorspace::lighten(map_pal_mtx[ , 3], .1)
map_pal_mtx[ , 2] <- colorspace::lighten(map_pal_mtx[ , 2], .2)
map_pal_mtx[ , 1] <- colorspace::lighten(map_pal_mtx[ , 1], .3)
map_pal_mtx[1, 1] <- '#ffffee'

map_pal <- as.vector(map_pal_mtx) %>% setNames(names(map_pal_raw))

panel_lgd_raw <- bi_legend(pal = map_pal, dim = 4,
                           x = 'Climate',
                           y = 'Non-climate')

si_fig <- panel_lgd_raw + 
  geom_text(data = cc_noncc_summary %>%
              filter(cc_spp + noncc_spp < 6), 
            aes(x = cc_spp, y = noncc_spp, label = lbl),
            hjust = .5, vjust = .5, color = 'black', size = 3) + 
  geom_text(data = cc_noncc_summary %>%
              filter(cc_spp + noncc_spp >= 6), 
            aes(x = cc_spp, y = noncc_spp, label = lbl),
            hjust = .5, vjust = .5, color = 'white', size = 3) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.title = element_text(hjust = 0, size = 12, face = 'bold', family = 'Arial'),
        panel.border = element_rect(color = 'black', fill = NA))
si_fig

ggsave(filename = 'figS3A_spp_legend_w_pcts.png', plot = si_fig, height = 3, width = 3, units = 'in', dpi = 300)

panel_d_biplot <- ggplot() +
  ### plot data:
  geom_raster(data = cc_noncc_spp_qtile_df, aes(x, y, fill = bi_class), show.legend = FALSE) +
  ### continents:
  geom_sf(data = land_sf,
          fill = 'grey96', color = NA, 
          size = .1) +
  geom_sf(data = globe_border,
          fill = NA, color = 'grey70', 
          size = .1) +
  # geom_rect(data = inset_bbox_df, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
  #           color = 'red', fill = NA, linewidth = 0.2) +
  bi_scale_fill(pal = map_pal, dim = 4) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_void() +
  theme(axis.text = element_blank(), 
        axis.title = element_blank())

# panel_d_biplot

3.2 Assemble figure

fig2_file <- here('5_ms_figs/fig2_mapping_impacts_spp.png')
panel_a_map  <- get_panel(panel_a_chi)
panel_b_map  <- get_panel(panel_b_cc)
panel_c_map  <- get_panel(panel_c_noncc)
panel_d_map  <- get_panel(panel_d_biplot)
panel_d_lgd <-  get_panel(panel_lgd_raw)
panel_abc_lgd <- get_legend(panel_a_chi)
panel_e_inset <- get_panel(inset_a_i)
panel_f_inset <- get_panel(inset_a_ii)
panel_g_inset <- get_panel(inset_a_iii) 

### parameterize panels for easier changing
fig_h <- 8.5
fig_w <- 6.25
ar    <- fig_h / fig_w
p_h <- 0.24
p_w <- p_h * 2 * ar
i_h <- p_h * 0.8
i_w <- i_h * ar

a <- c(x = 0.00, y = 0.75)
b <- c(x = 0.00, y = 0.50)
c <- c(x = 0.00, y = 0.25)
d <- c(x = 0.00, y = 0.0)
e <- c(x = .968 - i_w, y = 0.65)
f <- c(x = .968 - i_w, y = 0.45)
g <- c(x = .968 - i_w, y = 0.25)

fig2_panels <- ggdraw() +
  draw_plot(panel_a_map,  x = a['x'], y = a['y'], height = p_h, width = p_w) +
  draw_plot(panel_b_map,  x = b['x'], y = b['y'], height = p_h, width = p_w) +
  draw_plot(panel_c_map,  x = c['x'], y = c['y'], height = p_h, width = p_w) +
  draw_plot(panel_d_map,  x = d['x'], y = d['y'], height = p_h, width = p_w) +
  draw_plot(panel_e_inset,  x = e['x'], y = e['y'], height = i_h, width = i_w) +
  draw_plot(panel_f_inset,  x = f['x'], y = f['y'], height = i_h, width = i_w) +
  draw_plot(panel_g_inset,  x = g['x'], y = g['y'], height = i_h, width = i_w)

fig2_labels <- fig2_panels +
  draw_label('A.', x = a['x'], y = a['y'] + p_h - .05,  hjust = 0, vjust = 1, 
             size = 12, fontfamily = 'Arial', fontface = 'bold') +
  draw_label('B.', x = b['x'], y = b['y'] + p_h - .05,  hjust = 0, vjust = 1, 
             size = 12, fontfamily = 'Arial', fontface = 'bold') +
  draw_label('C.', x = c['x'], y = c['y'] + p_h - .05,  hjust = 0, vjust = 1, 
             size = 12, fontfamily = 'Arial', fontface = 'bold') +
  draw_label('D.', x = d['x'], y = d['y'] + p_h - .05,  hjust = 0, vjust = 1, 
             size = 12, fontfamily = 'Arial', fontface = 'bold') +
  draw_label('E.', x = e['x'] - .01, y = e['y'] + i_h,  hjust = 1, vjust = 1, 
             size = 12, fontfamily = 'Arial', fontface = 'bold') +
  draw_label('F.', x = f['x'] - .01, y = f['y'] + i_h,  hjust = 1, vjust = 1, 
             size = 12, fontfamily = 'Arial', fontface = 'bold') +
  draw_label('G.', x = g['x'] - .01, y = g['y'] + i_h,  hjust = 1, vjust = 1, 
             size = 12, fontfamily = 'Arial', fontface = 'bold')

### parameterize legends for easier changing of label positions:
l1 <- c(x = p_w + .02, y = 0.85, h = .10, w = .97 - p_w)
l2 <- c(x = p_w + .048, y = 0.06, h = .12, w = .12 * ar)

fig2_legends <- fig2_labels +
  ### draw legends:
  draw_plot(panel_abc_lgd, x = l1['x'], y = l1['y'],
            height = l1['h'], width = l1['w']) +
  draw_plot(panel_d_lgd, x = l2['x'], y = l2['y'],
            height = l2['h'], width = l2['w']) +
  ### Legend label for CHI:
  draw_label('Cumulative human impacts\n(A-C, E-G)',
             x = l1['x'] + l1['w']/2, y = l1['y'] + l1['h'], 
             hjust = 0.5, vjust = 0.5, angle = 0,
             size = 10, color = 'black', fontfamily = 'Arial', fontface = 'bold') +
  ### Legend labels for biscale:
  draw_label('Impact\nquartiles (D)',
             x = l2['x'] + .01, y = l2['y'] + l2['h'], 
             hjust = 0, vjust = 0, angle = 0,
             size = 10, color = 'black', fontfamily = 'Arial', fontface = 'bold') +
  draw_label('Climate →',
             x = l2['x'], y = l2['y'] - .005,
             hjust = 0, vjust = 1, angle = 0,
             size = 8, color = 'black', fontfamily = 'Arial', fontface = 'bold') +
  draw_label('Non-climate →',
             x = l2['x'] - .005, y = l2['y'],
             hjust = 0, vjust = 0, angle = 90,
             size = 8, color = 'black', fontfamily = 'Arial', fontface = 'bold')


ggsave(fig2_file, width = fig_w, height = fig_h, units = 'in', dpi = 300)
ggsave(fig2_file %>% str_replace('.png$', '.tif'), width = fig_w, height = fig_h, units = 'in', dpi = 300)

knitr::include_graphics(fig2_file)